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s_i ' Abstract 
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We introduce a modified algorithm to perform nonlinear filtering of a time series 
by locally linear phase space projections. Unlike previous implementations, the algo- 
rithm can be used not only for a posteriori processing but includes the possibility to 
perform real time filtering in a data stream. The data base that represents the phase 
space structure generated by the data is updated dynamically. This also allows filtering 
^ • of non- stationary signals and dynamic parameter adjustment. We discuss exemplary 

OO ' applications, including the real time extraction of the fetal electrocardiogram from ab- 

, dominal recordings. 
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Nonlinear projective filtering of time series is one of the most practical outcomes of the 
J^ 1 - theory of nonlinear dynamical systems as applied to real word observations. Signals with non- 

linear dynamical correlations are often not handled properly by spectral processing methods. 
Nonlinear dynamical systems exhibiting deterministic chaos have been proposed as an al- 
ternative paradigm for the study of such complex sequences. However, for the derivation 
of time series methods within this framework, heavy use had to be made of the theoretical 
properties of deterministically chaotic systems. This seems to severely restrict the range of 
systems they may be applied to - after all, very few real world phenomena are adequately 
described by a purely deterministic time evolution. Also nonlinear filtering procedures were 
designed to exploit the specific structure generated by deterministic dynamics in phase space. 
It turns out, however, that careful use of these algorithms can give excellent results also in 
situations where pure determinism is violated. The reason is that when represented in a 
low-dimensional phase space, also non-deterministic systems may exhibit structures suitable 
for filtering purposes. One example is the human electro-cardiogram (ECG), a signal that 
shows features of nonlinear determinism but also essential stochastic components. 

An introduction to deterministic chaos and dynamical systems as a framework for the 
analysis of time series recordings can be found in the monographs Refs. [Q, |2|. Review articles 
on nonlinear filtering algorithms containing further references to the original literature are 
Refs. [§ f§. How and why the ECG can be processed with nonlinear phase space filters is 
discussed in ||. 

Apart from these fundamental considerations, there are also some practical issues that 
have so far hampered widespread use of nonlinear filters. All of the methods that have been 
proposed in the literature are formulated as a posteriori filters. The whole signal has to be 
available before a cleaned version can be computed. The actual computation is invariably 
quite computer time intensive. One class of methods uses a global nonlinear function to 
represent an approximation to the dynamics. This function has to be determined by a 
delicate fitting procedure and the actual filtering scheme (for example Ref. ||) consists of an 



1 



iterative minimisation procedure in a high-dimensional space. The other class of algorithms 
approximates the dynamics or geometry in phase space by locally linear mappings. Here, 
covering phase space with small neighbourhoods is the most time consuming step, along 
with the need to solve a least squares or eigenvalue problem in each of these neighbourhoods. 
With fairly low-dimensional signals and small noise levels, fast neighbour search algorithms 
(see [Q for an overview) are very helpful in this regard. In any case, the posterior nature and 
the computational effort has been one of the major reasons why nonlinear phase space filters 
haven't seen more widespread use. The purpose of this article is to introduce modifications 
to a locally projective noise reduction scheme that make its use in real time signal processing 
feasible. 

In brief, we implement three main modifications to the locally projective noise reduction 
scheme that has been introduced in Ref. ||. Exactly the same strategy can be applied to 
modify the algorithms by Sauer || or that by Cawley and Hsu (To). (1) The data base of 
local neighbours that is needed to approximate the dynamics is restricted to points in the 
past, thereby rendering the filter causal. Of course, at the beginning no data base is available 
and noise reduction gradually becomes more effective as new points are collected. As a side 
effect, the curvature correction proposed in and discussed in Q can be carried out during 
the first sweep through the data. This will be explained below. (2) The number of neighbours 
required to carry out the correction for each point is limited to a number that is just sufficient 
for statistical stability of the local linear fits. (3) The last and most severe modification uses 
the fact that the dynamics is supposed to vary smoothly in phase space. Therefore, instead of 
determining the locally linearised dynamics at each point, the linear approximation is stored 
only for a collection of representative points which densely cover the observed set of points. 
Consequently, the local linear problem has to be solved only for points which are about to 
become a representative. 

The algorithm we present in this paper is based on the noise reduction scheme introduced 
in Ref. An alternative derivation can be found in Ref. 111]]. We refer the reader to 
these references as well as the monograph Ref. || for the motivation of the approach by the 
theory of deterministic dynamical systems. In a more general context, a scalar time series 
{s„},n = 1, . . . , N can be unfolded in a multi-dimensional effective phase space using time 
delay coordinates s„ = (s„_( m _iv, . . . , s n ) (r is a delay time). If {s n } is a scalar observation 
of a deterministic dynamical system, it can be shown under certain genericity conditions |l2| , 
[l3| that the reconstructed point set is a one-to-one image of the original attractor of the 
dynamical system. We will not explicitly assume here that there is such an underlying 
deterministic system. Nevertheless, general serial dependencies among the {s n } will cause 
the delay vectors {s„} to fill the available m-dimensional space in an inhomogeneous way. 
Linearly correlated Gaussian random variates will for example be distributed according to an 
anisotropic multivariate Gaussian distribution. Linear geometric filtering in phase space seeks 
to identify the principal directions of this distribution and project onto them. This concept is 
implemented by the singular systems approach, see Refs. jli], |l5|, |l6| and many others. The 
present algorithm can be seen as a nonlinear generalisation of this approach which takes into 
account that nonlinear signals will form curved structures in delay space. In particular, noisy 
deterministic signals form smeared-out lower-dimensional manifolds. Nonlinear phase space 
filtering seeks to identify such structures and project unto them in order to reduce noise. 

Thus, conceptually, the noise reduction algorithm consists of three main steps. (1) Find 
a low-dimensional approximation to the "attractor" described by the trajectory {s„}. (2) 
Project each point s„ on the trajectory orthogonally onto the approximation to the attractor 
to produce a cleaned vector s n (3) Convert the sequence of cleaned vectors s n back into the 
scalar time domain to produce a cleaned time series s n . 

The low-dimensional approximation to the point set can be constructed locally in delay 
coordinate space using a procedure that is very similar to a local version of principal com- 
ponent analysis. In order to define a neighbourhood in phase space around a point s n , find 
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Figure 1: Local linear approximation to a one-dimensional curve. Left: approximations are 
not tangents but secants and all the centres-of-mass (crosses) of different neighbourhoods are 
shifted inward with respect to the curvature. Right: a tangent approximation is obtained 
by shifting the centre-of-mass outward with respect to the curvature. The open square 
denotes the average of the centres-of-mass of adjacent neighbourhoods, the filled square is 
the corrected centre-of-mass. 

all of the points that are within a distance e of s„. A set of nearby points can be defined as 

W&> = K: n-An<n'<n, \\s' n - s n \\ < e} (1) 

where ||s^ — s„|| is the phase space distance between and s„. (We use the max norm to 
measure distances.) In this definition we already included two modifications with respect 
to ||. First, only points in the past (n' < n) are considered and second, neighbours are only 
considered that have been observed no longer than An time steps ago. The neighbourhood 
which is actually used for noise reduction is determined by finding the largest An < n 
for which the number of vectors in U^ n \ does not exceed a specified £/max- 

The local structure of the point set is approximated by a linear subspace formed essentially 
by local principal components. The most straightforward implementation is to compute the 
local centre-of-mass 

(8)W = |W(«)|- 1 Y, s «' ( 2 ) 

s„,€UW 

where ^ s , £W (») means summation over all |W^ n ^| vectors in U^ n >. Then the principal com- 
ponents can be calculated around this point. In that case, however, the linear subspace is 
not tangent to the curved manifold but rather intersects with it, as illustrated in Fig. |l|. 
Therefore, it is preferable to use a corrected centre-of-mass sr™-* given by 

S< n > =2(s)W _ \UW\- 1 Y ( s ) {n ' } ( 3 ) 

In the original implementation [||, neighbourhoods could contain future points. Thus, Eq.(^) 
could only be evaluated after one complete sweep through the data in which all the local 
centres-of-mass (s)( n ' were determined. In order to compute the tangent points s^ n \ a second 
sweep through the data set was necessary. In this modified implementation, the centre-of- 
mass vectors (s)'"' are stored when the point with index n is processed. Thus all these 
vectors are available for points with index n' < n. Thus, s^™-* can be formed immediately. 
Now, the local weighted covariance matrix 

C\f= [R(sn'-§ (n) )]i[R(s„'-I (n) )]i (4) 

is computed, where [•], denotes the i-th component of the vector in brackets. As discussed 
in Ref. |J, the weight matrix R is chosen to be diagonal with Rn and R mm large and 
all other diagonal entries Ru = 1. Now determine the orthonormal eigenvectors c q and 
eigenvalues of C|" using standard matrix techniques. A Q-dimensional manifold is then 
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locally approximated by those Q eigenvectors with the largest eigenvalues. The projected 
vector s„ is then given by: 

Q 

s n = + R 1 c? [ c? ' R ( s « " ■ ( 5 ) 

9=1 

In order to translate {s n } back into a scalar signal, we note that each scalar measure- 
ment s n appears as a component in m embedding vectors, s„, . . . , s„+ m _i. The corrected 
scalar time scries values s n are thus obtained by averaging the corresponding components of 

Sn , ■ ■ ■ , Sn+m-1 • 

For real time application, this means that the corrected value s n cannot be available 
before s n +m-i has been measured and processed. Usually, however, this delay window is a 
very short time, at least compared to the duration of the recording, and the procedure can 
be regarded as effectively on-line as long as the computations necessary to obtain s n can be 
carried out fast enough. 

So far, we have turned the procedure into a causal filter by restricting neighbour search 
to points defined by measurements made in the past. By further limiting the number of 
neighbours searched for, we have sped up the formation of the local covariance matrices 
considerably. Still, as the algorithm stands, we have to solve an (m x m) eigenvalue problem 
for each point that is to be processed. A large fraction of this work can be avoided on 
the base of the assumption that the local linear structure changes smoothly over phase 
space. By making local linear approximations we have already assumed smoothness of the 
underlying manifold in the C\ sense. In most physical systems, the additional assumption 
of C2 smoothness is not less justified. We cannot, however expect that the vectors which 
span the principal directions vary slowly from point to point. The reason is that often 
some eigenvalues are nearly degenerate and change indices from point to point when they 
are ordered by their magnitude. We thus refrain from interpolating principal components 
between phase space points. Instead, we choose a length scale h in phase space which is 
small enough such that the linear subspaces spanned by the local principal components can 
be regarded as effectively the same. Now we successively build up a data base of representative 
points for which the local points of tangency sS n ) , and the local principal directions c q , q = 
1,...,Q have been determined already. For each new point s„ that is to be processed, we go 
through this collection of points to determine whether a representative is available closer than 
h. In that case, we use the stored tangent point and principal directions of the representative 
in order to perform the projections. If not, a neighbourhood is formed around s n in which 
the eigenvalue problem is solved. The point s„ is then included in the list of representatives. 
By this procedure, all parts of the available space that is visited by observations are covered 
by representatives with a maximal distance h. If the dynamics and thus the geometry in 
phase space undergoes some change during the recording, new areas are visited which are 
automatically covered by new representatives. It is also possible to delete representatives 
which arc older than a given time span. A realisation of the algorithm that implements all 
the modifications discussed in this paper is publicly available as part of the TISEAN software 
project 0. 

As a first example, we show in Fig. ^| the result of applying the described procedure to a 
data set from an NMR laser experiment [pl| . The same data has also been used in Ref . pl| . 
The laser is periodically driven and once per driving cycle the envelope of the laser output 
is recorded. The resulting sampling rate is 91 Hz. At this rate, the modified nonlinear noise 
reduction scheme can be easily carried out in real time on a Pentium II processor at 200 MHz. 
Since further iterations can be carried out after the time corresponding to one embedding 
window without interfering with the previous steps, we could perform up to three iterations 
on a dual Pentium II workstation at 300 MHz in real time. The figure shows the result after 
two iterations. Projections from m = 7 down to Q = 2 dimensions were used, at least 100 
neighbours were requested at e = 200 A/D units. The history was limited to 20000 samples, 
or 220 s. This fairly large data base is needed since the initial noise level is already small 
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Figure 2: Result of nonlinear filtering of an NMR laser time series. An enlargement of about 
one quarter of the total linear extent of the attractor is shown. 
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Figure 3: Result of nonlinear filtering of an MCG time series. Note that the noise is not white. 
It contains for example contributions from non-cardiac muscle activity and fluctuations in 
the magnetic background field. 

(less than 2% [^9|) and small neighbourhoods are required to avoid curvature artefacts. The 
maximal distance of representative points was chosen to be h = 120 A/D units. 

The actual acceleration resulting from the above modifications strongly depends on the 
situation and it is difficult to give general rules and benchmarks. Let us however study a 
realistic example in some detail to illustrate the main points. In electrophysiological research, 
it is quite attractive to augment the measurements of electric potentials with recordings of 
the magnetic field strength. The latter penetrate intervening tissue much more efficiently. 
Of particular interest are magneto- encephalographic (MEG) recordings which allow to access 
regions of the brain noninvasively which cannot be monitored electrically using surface elec- 
trodes. In cardiology, the magneto- cardiogram (MCG) provides additional information to the 
traditional electrocardiogram (ECG). A particular application is the noninvasive monitoring 
of the fetal heart which is otherwise complicated by shielding of the electric field by inter- 
vening tissue. A common problem with magnetic recordings, however, is that the fields are 
rather feeble and the measurements have to be carried out in a shielded room. Even then, 
noise remains a major challenge for this experimental technique. We will demonstrate in the 
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Figure 4: Result of nonlinear filtering of an MCG time series. Same data as Fig. || but in 
time delay representation (enlargement). 



method 



CPU time 



a) 
b) 
c) 
d) 
e) 
f) 



all neighbours 

box assisted neighbour search 
all neighbours in past 
n — n! < 5 s 

(d) and [/max < 200 

(e) and reuse of representatives 



165 s 
31 s 
82 s 
64 s 
15 s 
2 s 



Table 1: Computation time for one iteration of the nonlinear noise reduction scheme, applied 
to 10 s of an MCG recording sampled at 1000 Hz. See text for details. 

following how nonlinear noise reduction could be used for continuous MCG monitoring. 

We use an MCG recording of a normal human subject at rest. The data was kindly 
provided by Carsten Sternickel at the University of Bonn. The sampling rate was 1000 Hz, 
which is quite high. For the signal processing task, any sampling rate above about 200 Hz 
would be sufficient. (Below 200 Hz, the spike representing the depolarisation of the ventricle 
might not be resolved properly). In order to cover a significant fraction of one cardiac 
cycle by an embedding window, we chose an embedding with delay r = 10 ms in m = 10 
dimensions. Neighbourhoods were formed with a radius of e = 0.1 (in the uncalibrated A/D 
units of the recording), about three times the noise level estimated by visual inspection. In 
Figs. H and |], the result of two iterations of the noise reduction scheme is shown. We quote 
computation times for the processing of 10 s of MCG on a Pentium II processor at 300 MHz, 
determined on a PC running the Linux operating system. The timing results are summarised 
in Table [|. The data base of representatives for (f) was formed by assuming that the local 
linear subspaces are equivalent on length scales of the order of 0.06 A/D units. 

Like in the case of the MCG time series demonstrated above, nonlinear phase space filter- 
ing can be successfully applied to ECG recordings pp| . The modifications of the algorithm 
described in the present paper allow to perform the filtering as a real time application on 
common personal computers without a noticeable impairment in signal quality compared to 
previous implementations. Since this application is similar to the MCG filtering described 
above, we do not need to state further details. Instead, let us discuss a related signal sep- 
aration problem of practical relevance, the extraction of the fetal ECG from abdominal 
recordings during pregnancy [^l], ^2| . In order to obtain the strongest possible signal of the 
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Figure 5: Fetal ECG extraction by nonlinear filtering. Upper trace: the original abdominal 
recording. Lower trace: the extracted fetal ECG. See text for details. 

fetal cardiac activity, electrodes are usually placed on the maternal abdomen. The upper 
trace of Fig. |^ contains such an abdominal recording. (In this particular case one abdominal 
and one cervical electrode has been used; the data were kindly provided by J. F. Hofmeister, 
Denver |^|.) Due to the large size of the maternal heart in comparison to the fetal one, the 
maternal ECG represents the dominant signal component. Despite a strong noise component, 
the spike-like ventricular complexes of the fetal signal are present. The noise floor is due to 
action potentials of the muscular tissue surrounding the electrodes and other measurement 
noise. (Note the small overall amplitude of the signal.) The nonlinear filter described in this 
paper can be used to extract the fetal signal in a two-step procedure: in the first step the 
maternal signal is cleaned up by considering the noise level to be of the magnitude of the 
fetal spikes. This identifies the noise and the fetal component together as a contamination 
of the maternal waveform. The difference between the abdominal recording and the clean 
maternal component then yields a noisy fetal signal. The fetal ECG can be separated from 
the noise in a second nonlinear noise reduction step. The result of this procedure is shown 
in Fig. H The abdominal recording (upper trace) is processed in two steps to yield a fairly 
clean fetal ECG (lower trace). The algorithm described in this paper allows the extraction 
of the fetal ECG as a real time application on a Laptop PC (Pentium processor at 133 MHz, 
Linux operating system) for a sampling rate of 250 Hz. 

In conclusion, we have demonstrated that by certain modifications to the algorithms de- 
scribed in the literature, nonlinear projective noise reduction can be turned into a signal 
processing tool that can in many situations run in real time in a data stream. Necessary 
ingredients are (1) the formulation of the algorithm as a causal filter, relying only on informa- 
tion that is available at recording time, (2) a general speed-up of the procedure by restricting 
neighbour search to the immediate past (at the expense of peak noise reduction performance), 
and (3) further speed-up by using a data base of previously processed phase space points. By 
putting a time restriction on the data base for the formation of local subspaces, processing 
of non-stationary signals with slowly drifting parameters becomes also possible. 

Carsten Sternickel was so kind to let us use his MCG recordings. Leci Flepp provided 
the NMR laser time series and John F. Hofmeister made his fetal ECG recordings available. 
We thank Daniel Kaplan, Rainer Hegger, and Holger Kantz for useful discussions. This work 
was supported by the SFB 237 of the Deutsche Forschungsgemeinschaft. 
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